Sujet2-Etape1

ℹ Using "','" as decimal and "'.'" as grouping mark. Use `read_delim()` for more control.

1 Etape 1 : Première manipulation du RPG

1.1 Je choisis un point sur la carte : https://www.google.fr/maps

=> clic droit sur la carte + clic gauche sur les coordonnées (copiées dans le presse-papier)

J’indique un rayon en mètres (< 15000) pour sélectionner les parcelles autour du point

Code

# coord_gmaps<- rstudioapi::showPrompt(title = "Collez les coordonnées", message = "coordonnées Gmaps", default = "")
coord_gmaps<-"46.46946179131805, -1.4932577154775046"
lat<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,1])
lon<-as.numeric(str_split(coord_gmaps,fixed(","),simplify = TRUE)[,2])
# rayon<-as.numeric(rstudioapi::showPrompt(title = "Rayon", message = "Rayon (en m)", default = ""))
rayon<-5000

# création d'une table sf «point» avec les coordonnées saisies
# transformation des coordonnées en syst de proj 2154 (Lambert II - Français) 
point<-data.frame(lon,lat,rayon) %>% 
  st_as_sf(coords = c("lon","lat"),crs = "EPSG:4326") %>%
  mutate(coord_pt_gps=st_as_text(geometry)) %>% 
  st_transform("EPSG:2154") %>% 
  st_sf() %>% clean_names() %>% 
  rename(geom=geometry)

# st_crs(point)
# st_geometry_type(point)
# plot(point)

La table des parcelles agricoles 2021 se trouve sur un serveur PostgreSQL/PostGIS.

1.2 Je me connecte au serveur PostgreSQL avec le mot de passe disponible sur Onyxia/Mes services/PostgreSQL/Read me

Code

# le mot de passe est stocké dans un secret Vault

# postgresql_password <- rstudioapi::askForPassword(prompt = "Entrez le password PostgreSQL")
postgresql_password <- "1tfawt3nj7fgzo3w7cma"
  
# Connection à PostgreSQL
cnx <- dbConnect(PostgreSQL(),
                user = "projet-funathon",
                password = postgresql_password,
                host = "postgresql-438832",
                dbname = "defaultdb",
                port = 5432,
                options="-c search_path=rpg,public") # specify what schema to connect to

1.3 je lance une requête SQL pour sélectionner les parcelles situées autour de mon point dans le rayon choisi.

Code

# suppression de la table «point» si elle existe
dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.point CASCADE;")

<PostgreSQLResult>

Code

# écriture de la table point dans la base PostGis
write_sf(point, cnx, append = T)

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.point ADD CONSTRAINT point_pkey PRIMARY KEY(coord_pt_gps);")

<PostgreSQLResult>

Code

# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

<PostgreSQLResult>

Code

# dbExecute(cnx,"CREATE INDEX ON rpg.point USING gist(geom);")

# 3 - Exécution de la requête de découpage du RPG autour du point sur PostGis  -------

dbSendQuery(cnx,"DROP TABLE IF EXISTS rpg.parc_prox CASCADE;")

<PostgreSQLResult>

Code

dbSendQuery(cnx,"CREATE TABLE rpg.parc_prox AS
    SELECT row_number() OVER () AS row_id, p.coord_pt_gps, p.rayon, r.*  
    FROM rpg.point p, rpg.parcelles r 
    WHERE ST_DWithin(p.geom,r.geom,p.rayon);")

<PostgreSQLResult>

Code

# ajout d'une clé primaire
dbSendQuery(cnx,"ALTER TABLE rpg.parc_prox ADD CONSTRAINT parc_prox_pk PRIMARY KEY(id_parcel);")

<PostgreSQLResult>

Code

# ajout d'un index
dbSendQuery(cnx,"CREATE INDEX parc_prox_geom_idx ON rpg.parc_prox USING gist(geom);")

<PostgreSQLResult>

Code

# 4 - Téléchargement des parcelles proches sous R-------------------------------

# parc_prox<-read_sf(cnx,parc_prox)

parc_prox<-st_read(cnx, query="select * from rpg.parc_prox;")

1.4 J’affiche les parcelles autour de mon point avec une carte interactive leaflet

Code

#plot(st_geometry(parc_prox))
#plot(st_geometry(point),add = T,col = "red")

# Transformation de la projection car leaflet ne connait que le WGS 84
carte_parc_prox<- parc_prox %>% st_transform(4326)

# Marqueur du point
pt_mark<- point %>% st_transform(4326)

# ajout du libellé des cultures
carte_parc_prox_lib<-carte_parc_prox %>% 
  left_join(lib_cult %>% select(-code_groupe_culture),by=c("code_cultu"="code_culture")) 
#création d'un label ad hoc à afficher en surbrillance au passage de la souris sur la carte.
labels<-sprintf("<strong>id_parcel : </strong>%s<br/>
                <strong>Groupe culture : </strong>%s<br/>
                <strong>Culture : </strong>%s<br/>
                <strong>Surface (ha) : </strong>%s<br/>
                <strong>Département : </strong>%s<br/>
                <strong>Commune : </strong>%s<br/>",
                parc_prox$id_parcel,carte_parc_prox_lib$libelle_groupe_culture,
                carte_parc_prox_lib$libelle_culture,parc_prox$surf_parc,
                parc_prox$insee_dep,parc_prox$nom_com) %>% 
  lapply(htmltools::HTML)

# labels

# création d'une palette de couleurs associée au groupe de culture
factpal <- colorFactor("Paired", parc_prox$code_group)

carte_parc_prox_html <- leaflet(carte_parc_prox_lib) %>% 
   # addProviderTiles("Esri.WorldImagery") %>%
   # addTiles() %>% 
  addTiles("http://wxs.ign.fr/essentiels/wmts?REQUEST=GetTile&SERVICE=WMTS&VERSION=1.0.0&STYLE=normal&TILEMATRIXSET=PM&FORMAT=image/jpeg&LAYER=ORTHOIMAGERY.ORTHOPHOTOS&TILEMATRIX={z}&TILEROW={y}&TILECOL={x}") %>%
  addPolygons( #fillColor="white",
               fillColor=~factpal(code_group),
               weight = 2,
               opacity = 1,
               color = "#ffd966",
               dashArray = "3",
               fillOpacity = 0.5,
               highlight = highlightOptions(
                 weight = 5,
                 color = "#A40000",
                 dashArray = "",
                 fillOpacity = 0.0,
                 bringToFront = TRUE),
               label = labels,
               labelOptions = labelOptions(
                 style = list("font-weight" = "normal", padding = "3px 8px"),
                 textsize = "15px",
                 direction = "auto",
                 encoding="UTF-8")) %>% 
  addMarkers(data=pt_mark,~lon,~lat, popup = ~coord_pt_gps, label = ~coord_pt_gps)

Warning in RColorBrewer::brewer.pal(max(3, n), palette): n too large, allowed maximum for palette Paired is 12
Returning the palette you asked for with that many colors

Code

# leaflet(pt_mark) %>% addTiles() %>% addMarkers(~lon,~lat)

carte_parc_prox_html

Code

# saveWidget(widget = carte_parc_prox_html, file = "carte_parc_prox_html.html")

1.5 Quelle est la structure des parcelles agricoles autour de mon point ?

Code

t1 <- parc_prox %>% st_drop_geometry() %>% count(code_group) %>% 
  add_tally(n) %>% 
  mutate(n_pct=round(100*n/nn,1)) %>% 
  select(-nn) %>% rename(n_parcelles=n) %>%
  # adorn_totals() %>% 
cbind(
  # comptage des surfaces
parc_prox %>% st_drop_geometry() %>% count(code_group,wt=surf_parc) %>% 
      add_tally(n) %>% 
      mutate(surf_pct=round(100*n/nn,1)) %>%
      select(-nn) %>%  
      rename(surf_parc_ha=n) %>% select(surf_parc_ha, surf_pct) # %>% 
       # adorn_totals()
) %>% left_join(lib_group_cult,by=c("code_group"="code_groupe_culture")) %>% 
  select(code_group,libelle_groupe_culture,everything()) %>% 
  # arrange(as.numeric(code_group)) %>% 
  arrange(desc(surf_parc_ha)) %>% 
  adorn_totals() %>% 
  mutate(taille_moy_parc=round(surf_parc_ha/n_parcelles,1))

Storing counts in `nn`, as `n` already present in input
Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.

Code

t1  %>% 
  setNames(c("code","groupe de cultures","nombre de parcelles","(%)","surface (ha)","surface (%)","taille moyenne (ha)")) %>% 
  kable(
    format="html",
    caption="<span style='font-size:medium'>Groupes de cultures <strong>locales</strong> par surfaces décroissantes</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(t1), bold = T, color = "white", background = "grey")

Structure des cultures au niveau local

code	groupe de cultures	nombre de parcelles	(%)	surface (ha)	surface (%)	taille moyenne (ha)
2	Maïs grain et ensilage	306	18,8	1 669,06	30,4	5,5
18	Prairies permanentes	402	24,7	1 168,78	21,3	2,9
1	Blé tendre	178	10,9	737,43	13,4	4,1
19	Prairies temporaires	175	10,7	501,54	9,1	2,9
4	Autres céréales	88	5,4	428,99	7,8	4,9
6	Tournesol	76	4,7	319,70	5,8	4,2
16	Fourrage	39	2,4	141,23	2,6	3,6
5	Colza	31	1,9	133,49	2,4	4,3
7	Autres oléagineux	13	0,8	74,82	1,4	5,8
11	Gel (surfaces gelées sans production)	86	5,3	73,82	1,3	0,9
25	Légumes ou fleurs	22	1,4	70,15	1,3	3,2
28	Divers	180	11,1	65,18	1,2	0,4
3	Orge	20	1,2	52,31	1,0	2,6
15	Légumineuses à grains	9	0,6	30,65	0,6	3,4
8	Protéagineux	3	0,2	19,63	0,4	6,5
Total	-	1 628	100,1	5 486,78	100,0	3,4

Code

# rm(t1)

1.6 Comparaison avec la répartition des cultures au niveau départemental et national

Code

# jointure spatiale avec la couche département pour récupérer le département où tombe le point

# ouvrir la couche des communes (à convertir en Lambert II epsg 2154)
dep <- s3read_using(
    FUN = sf::read_sf,
    layer = "departement",
    object = "2023/sujet2/diffusion/ign/adminexpress_cog_simpl_000_2023.gpkg",
    bucket = "projet-funathon",
    opts = list("region" = "")) %>% 
  st_transform(2154)

# jointure
df<-point %>% st_join(dep) %>% st_drop_geometry() %>% select(insee_dep)
dep_pt<-df[1,1]
rm(df)

# sinon sélection du département regroupant la plus grande surface agricole
# cas où le cercle recouvre plusieurs départements
# df<-parc_prox %>% st_drop_geometry() %>% 
#   count(insee_dep,wt=surf_parc) %>% 
#   arrange(desc(n))
# dep_pt<-df[1,1]
# rm(df)

# calcul des % surfaces autour du point
stat_pt <- parc_prox %>% st_drop_geometry() %>% 
  count(code_group,wt=surf_parc) %>% add_tally(n) %>% 
  mutate(pct_surf_local=round(100*n/nn,1)) %>%
  select(code_group, pct_surf_local) 

Storing counts in `nn`, as `n` already present in input
ℹ Use `name = "new_name"` to pick a new name.

Code

# récup des % surfaces départementales
stat_dep_pt<-s3read_using(
  FUN=readr::read_rds,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_by_dep.rds",
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
  filter(insee_dep %in% dep_pt) %>% 
  select(insee_dep,code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_dep = pct_surf)

# récup des % surfaces nationales
stat_fm<-s3read_using(
  FUN=readr::read_csv,
  object = "2023/sujet2/diffusion/resultats/stat_group_cult_fm.csv",
  col_types = cols(code_group = col_character()),
  bucket = "projet-funathon",
  opts = list("region" = "")) %>% 
 select(code_group,libelle_groupe_culture,pct_surf) %>% 
  rename(pct_surf_fm = pct_surf)

# appariement des stas locales, départementales, nationales
stat_compar<-stat_fm %>% 
  left_join(stat_dep_pt %>% select(code_group, pct_surf_dep),by="code_group") %>% 
  left_join(stat_pt,by="code_group") %>% 
  select(libelle_groupe_culture,pct_surf_local,pct_surf_dep,pct_surf_fm) %>% 
  arrange(desc(pct_surf_local)) %>% adorn_totals() 

stat_compar %>% 
  setNames(c("Groupe de cultures","surf. locales (%)","surf. départ. (%)","surf. France m. (%)")) %>% kable(
    format="html",
    caption="<span style='font-size:medium'>Comparaison des surfaces locales, départemenales et nationales</span>",
    format.args = list(decimal.mark = ",", big.mark = " "),
    booktabs = TRUE) %>%
  kable_styling(font_size = 15) %>% 
  gsub("font-size: initial !important;",
       "font-size: 20pt !important;",.)%>%
  kable_styling(bootstrap_options = c("striped", "hover")) %>% 
  row_spec(nrow(stat_compar), bold = T, color = "white", background = "grey")

Structure des cultures

Groupe de cultures	surf. locales (%)	surf. départ. (%)	surf. France m. (%)
Maïs grain et ensilage	30,4	19,4	10,0
Prairies permanentes	21,3	29,0	27,7
Blé tendre	13,4	16,4	17,7
Prairies temporaires	9,1	9,5	5,3
Autres céréales	7,8	8,1	3,9
Tournesol	5,8	3,7	2,5
Fourrage	2,6	4,1	3,7
Colza	2,4	2,3	3,5
Autres oléagineux	1,4	0,4	0,7
Gel (surfaces gelées sans production)	1,3	0,6	1,6
Légumes ou fleurs	1,3	1,2	1,5
Divers	1,2	1,0	1,3
Orge	1,0	2,3	6,2
Légumineuses à grains	0,6	0,4	0,2
Protéagineux	0,4	0,9	1,2
Plantes à fibres	NA	0,3	0,4
Riz	NA	NA	0,0
Estives et landes	NA	0,1	8,1
Vergers	NA	0,1	0,4
Vignes	NA	0,2	2,2
Fruits à coque	NA	0,0	0,2
Oliviers	NA	NA	0,0
Autres cultures industrielles	NA	0,1	1,8
Total	100,0	100,1	100,1

1.7 Graphique de comparaison des cultures au niveau local, départemental et national

Code

# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# levels(as.factor(tab_piv$secteur))
# tab_piv$secteur<-factor(tab_piv$secteur, levels=c("france","departement","local"))
# tab_piv<-tab_piv %>% arrange(desc(secteur), desc(value))

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
p<-ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value, 
                       fill=factor(secteur, levels=c("france","departement","local")))) + 
  geom_col(position = "dodge") +
  labs(title="Surfaces comparées des 10 principales cultures locales, en %", x="Culture", y = "%", fill = "Secteur") +
  theme_classic()

# je flippe le graphique pour avoir des barres horizontales  
p+coord_flip()

différences de structure locale, departementale, nationale

1.8 Je teste un graphique par secteur avec facet_wrap

Code

# je sélectionne les 10 groupes de cultures les plus répandus au niveau local 
tab<-stat_compar %>% filter(libelle_groupe_culture!="Total") %>% slice_head(n=10) %>% 
  rename(local=pct_surf_local, departement=pct_surf_dep, france=pct_surf_fm)

# je transpose la table pour rassembler toutes les valeurs dans une seule variable value (ggplot oblige)
tab_piv<-tab %>% pivot_longer(!libelle_groupe_culture) %>% rename(secteur=name) 

# je mets à 0 les valeurs manquantes
#tab_piv<-tab_piv %>% coalesce(value,0L)  
tab_piv[is.na(tab_piv)] <- 0

# ggplot => problème de tri : affichage des prairies avant le maïs ?    
ggplot(tab_piv, aes(x = reorder(libelle_groupe_culture,+value), 
                       y = value)) + 
  geom_col(fill = "lightblue", colour = "black",position = "dodge") +
  labs(title="Surface par culture", x="Culture", y = "%", fill = "Secteur") +
  geom_text(aes(label = value), hjust = -0.3, size = 8/.pt, colour = "black") +
  theme_classic() + coord_flip() + 
  facet_wrap(~secteur,nrow=3,scales='free')

différences de structure locale, departementale, nationale

1.9 Pour aller plus loin, comparer la taille moyenne des parcelles selon les secteurs